1. Wczytanie zestawu danych punktowych oraz nadanie im ukladu wspolrzednych ETRS 1989 Poland CS2000 Zone 7 - EPSG:2178, nastepnie zapisanie ich do nowego pliku csv - zestaw1_epsg_2178.csv.
library(sp)

data_points <- read.csv('zestaw1.csv', colClasses = c('numeric', 'numeric'))
head(data_points)
##       Long      Lat
## 1 19.89514 50.07248
## 2 20.03497 50.07282
## 3 19.93090 50.04349
## 4 19.86676 50.06638
## 5 19.93647 50.05540
## 6 19.93338 50.06552
dim(data_points)
## [1] 2000    2
coord <- SpatialPoints(cbind(data_points$Long, data_points$Lat), 
                       proj4string = CRS('+proj=longlat'))
class(coord)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
length(coord)
## [1] 2000
head(coord)
## SpatialPoints:
##      coords.x1 coords.x2
## [1,]  19.89514  50.07248
## [2,]  20.03497  50.07282
## [3,]  19.93090  50.04349
## [4,]  19.86676  50.06638
## [5,]  19.93647  50.05540
## [6,]  19.93338  50.06552
## Coordinate Reference System (CRS) arguments: +proj=longlat +datum=WGS84
## +no_defs
coordUTM <- spTransform(coord, CRS('+init=epsg:2178'))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum European_Terrestrial_Reference_System_1989 in
## Proj4 definition
head(coordUTM)
## SpatialPoints:
##      coords.x1 coords.x2
## [1,]   7420912   5549067
## [2,]   7430922   5548966
## [3,]   7423426   5545805
## [4,]   7418870   5548418
## [5,]   7423844   5547124
## [6,]   7423639   5548253
## Coordinate Reference System (CRS) arguments: +proj=tmerc +lat_0=0
## +lon_0=21 +k=0.999923 +x_0=7500000 +y_0=0 +ellps=GRS80 +units=m
## +no_defs
length(coordUTM)
## [1] 2000
#uklad CRS pomyslnie zmieniony z WGS84 na EPSG:2178

coordUTM_df <- as.data.frame(coordUTM) #col.names nie dziala - zbadac, z jakiego powodu
colnames(coordUTM_df) <- c('Lon', 'Lat')
head(coordUTM_df)
##       Lon     Lat
## 1 7420912 5549067
## 2 7430922 5548966
## 3 7423426 5545805
## 4 7418870 5548418
## 5 7423844 5547124
## 6 7423639 5548253
dim(coordUTM_df) #dlugosc zgadza sie z dlugoscia wczytanych danych
## [1] 2000    2
write.csv(coordUTM_df, "zestaw1_epsg_2178.csv", row.names = FALSE)
2. Wczytanie danych w ukladzie wspolrzednych EPSG:2178 oraz pliku shapefile zawierajacego osiedla.
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(ggplot2)
library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
data <- read.csv('zestaw1_epsg_2178.csv')
head(data)
##       Lon     Lat
## 1 7420912 5549067
## 2 7430922 5548966
## 3 7423426 5545805
## 4 7418870 5548418
## 5 7423844 5547124
## 6 7423639 5548253
dim(data)
## [1] 2000    2
str(data) #2 zmienne - Lon oraz Lat numerical
## 'data.frame':    2000 obs. of  2 variables:
##  $ Lon: num  7420912 7430922 7423426 7418870 7423844 ...
##  $ Lat: num  5549067 5548966 5545805 5548418 5547124 ...
sum(is.na(data)) #brak pustych danych
## [1] 0
qplot(data$Lon, geom='histogram')
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

districts <- shapefile('osiedla.shp')
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in Proj4 definition: +proj=tmerc +lat_0=0 +lon_0=21 +k=0.999923 +x_0=7500000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs,
##  but +towgs84= values preserved
Wczytane dane maja rozmiar 2000 x 2, zaiweraja dlugosc i szerokosc geograficzna - typ numerical, zapisane sa w data.frame o nazwie data, brak pustych danych.

3. Prezentacja mapy dzielnic/osiedl w Krakowie z nalozonymi na nia wczytanymi punktami.
map_cracow <- ggplot() + 
  geom_polygon(data = districts, aes(x = long, y = lat, group = group), show.legend = FALSE, color = "white", fill = "darkgrey")  + coord_fixed()
## Regions defined for each Polygons
summary(districts)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##       min     max
## x 7413437 7443955
## y 5537344 5555031
## Is projected: TRUE 
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=21 +k=0.999923 +x_0=7500000 +y_0=0
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Data attributes:
##    GESTOSC_ZA         ID             IL_K_MOBIL         IL_K_NIEMO       
##  Min.   :    0   Length:141         Length:141         Length:141        
##  1st Qu.:  293   Class :character   Class :character   Class :character  
##  Median : 1625   Mode  :character   Mode  :character   Mode  :character  
##  Mean   : 5671                                                           
##  3rd Qu.: 7722                                                           
##  Max.   :30031                                                           
##   IL_K_POPRO         IL_K_PROD          IL_K_PRZED         IL_L_MOBIL       
##  Length:141         Length:141         Length:141         Length:141        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   IL_L_NIEMO         IL_L_POPRO         IL_L_PROD          IL_L_PRZED       
##  Length:141         Length:141         Length:141         Length:141        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   IL_M_MOBIL         IL_M_NIEMO         IL_M_POPRO         IL_M_PROD        
##  Length:141         Length:141         Length:141         Length:141        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   IL_M_PRZED          ILOSC_K_10     ILOSC_KOBI         ILOSC_MEZC       
##  Length:141         Min.   :  0.0   Length:141         Length:141        
##  Class :character   1st Qu.:103.3   Class :character   Class :character  
##  Mode  :character   Median :110.8   Mode  :character   Mode  :character  
##                     Mean   :113.2                                        
##                     3rd Qu.:119.6                                        
##                     Max.   :488.7                                        
##   ILOSC_MIES            KOD               MAPID              MSLINK         
##  Length:141         Length:141         Length:141         Length:141        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##   NAZWA_JEDN         NR_JEDN_UR         OZN_JEDN_U          POWIERZCHN      
##  Length:141         Length:141         Length:141         Min.   :   83366  
##  Class :character   Class :character   Class :character   1st Qu.:  892431  
##  Mode  :character   Mode  :character   Mode  :character   Median : 1638151  
##                                                           Mean   : 2318933  
##                                                           3rd Qu.: 3018523  
##                                                           Max.   :11796876  
##      WSP_X           WSP_Y        
##  Min.   :21980   Min.   :-301839  
##  1st Qu.:28779   1st Qu.:-295035  
##  Median :31681   Median :-292382  
##  Mean   :31075   Mean   :-291027  
##  3rd Qu.:33641   3rd Qu.:-286706  
##  Max.   :37118   Max.   :-275075
map_cracow + geom_point(data=data, aes(x=Lon, y=Lat), alpha=0.8, size = 1.2, colour="darkred") +ggtitle("Wykroczenia na mapie Krakowa\n z podzialem na dzielnice") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

4. Wydzielenie klastrow o zwiekszonej intensywnosci wystepowania wykroczen przy uzyciu algorytmow grupowania gestosciowego DBSCAN/ HBDSCAN / OPTICS.

4a
• DBSCAN (density-based spatial clustering of applications with noise) - nieliniowy algorytm uczenia sie bez nadzoru, ktory wykorzystuje technikę osiągalności i łączności gęstości. DBSCAN dzieli dane na grupy punktow majace wspolne cechy charakterystyki lub skupienia.

• Zasada działania algorytmu:
  1) Losowo wybierany jest punkt p
  2) Pobierane sa wszystkie punkty, ktorch gestosc jest osiagalna wzgledem punktu p na podstawue maksymalnego promienia sasiedztwa (EPS) oraz minimalnej liczby punktow w sasiedzwie EPS (MinPts)
  3) Dla kluczowych puntkow p tworzony zostaje klaster, w innym przypadku punkt p zostaje zaklasyfikowany jako punkt odstajacy lub szum.
  
• Klasyfikacja punktów:
  1) Centralny (core point) - przynajmniej MinPts liczba punktow (wlacznie z samym punktem p) znajduje sie w otoczeniu p z promieniem rownym EPS
  2) Graniczny (border point) - jezeli punkt jest osiagalny z punktu centralnego i w otoczeniu znajduje sie liczba punktow mniejsza, niz MinPts
  3) Odstajacy - jezeli punkt nie jest ani granicznym ani centralnym zotaje uznany za odstajacy
  
• Plusy DBSCAN:
  - dobre dzialanie przy dowolnych ksztaltach klastrow
  - odporny na wartosci odstajace, jest w stanie je wykryc
  - nie trzeba wczesniej z gory okreslac liczby klastrow
  
• Minusy DBSCAN:
  - charakterystyka tworzonych klastrow jest okreslana przez parametry MinPts oraz EPS przez co jezeli ciezko jest czasami utworzyc klastry o znaczacych roznicach w gestosci
  - czasami wyznaczenie odpowiedniej odleglosci sasiedztwa EPS nie nalezy do latwych i wymaga dodatkowej wiedzy
  
• Zastosowanie DBSCAN z parametrami odpowiednio MinPts i EPS rownymi:
  5/10, 10/50, 10/100, 10/200
library(dbscan)

dbscan_clusters1 <- dbscan(data, minPts = 5, eps = 10)
dbscan_clusters1
## DBSCAN clustering for 2000 objects.
## Parameters: eps = 10, minPts = 5
## The clustering contains 50 cluster(s) and 1569 noise points.
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 1569    7    8   10   13    5    7   11    5    8   10   11    6   14    9   12 
##   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31 
##   11   11   14   13   18    5   20    7    9    7    5   10    8    7    5    9 
##   32   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47 
##    6    6    5    8    6    6    6   13    5   16    7   11    9    5    6    5 
##   48   49   50 
##    5    6    5 
## 
## Available fields: cluster, eps, minPts
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = dbscan_clusters1$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=5, eps=10") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

dbscan_clusters2 <- dbscan(data, minPts = 10, eps = 50)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = dbscan_clusters2$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=10, eps=50") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

dbscan_clusters3 <- dbscan(data, minPts = 10, eps = 100)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = dbscan_clusters3$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=10, eps=100") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

dbscan_clusters4 <- dbscan(data, minPts = 10, eps = 200)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color =dbscan_clusters4$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("DBSCAN, minPts=10, eps=200") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

• Wnioski z DBSCAN:
1) Dla minPts = 5 i EPS = 10 bardzo mala ilosc klastrow, najwieksze skupisko w rejonie dzielnicy Stare Miasto, nastepnie mozemy zauwazyc klastry o zdecydowanie mniejszej gestosci w dzielniach Pradnik Bialy, Pradnik Czerwony, Mistrzejowice oraz Bienczyce. Szczatkowe ilosci punktow wystepuja rowniez na poludniu w dzielnicy Swoszowice.

2) Dla minPts = 10 i EPS = 50 zaczynamy zauwazac wyrazne klastry w okolicy dzielnicy Starego Miasta oraz pojedyncze wystepujace na terenie Nowej Huty/Wzgorz Krzeslawickich.

3) Dla minPts = 10 i EPS = 100 coraz wiecej klastrow o duzej gestosci w dzielnicy Starego Miasta, Grzegorzki, Czyzyny oraz Nowej Huty.

4) Dla minPts = 10 i EPS = 200 najwieksze, bardzo wyrazne skupisko klastrow w dzielnicy Stare Miasto, nastepnie mozemy zauwazyc wieksze klastry punktow w dzielnicy Krowodrza, Pradnik Bialy, Pradnik Czerwony, Mistrzejowice, Grzegorzki. Mniejsze klastry znajdziemy w Czyzynach, Wzgorzach Krzeslawickich oraz Nowej Hucie. Pojedynczy klaster wystepuje w Podgorzu Duchackim.


4b
• HDBSCAN (hierarchical density-based spatial clustering) - bardziej zaawansowana wersja algorytmu DBSCAN. Algorytm wykorzystuje podejscie oparte na gestosci - zamiast szukac klastrow posiadajacych okreslony ksztalt szuka regionow danych, ktore sa gestsze od otaczajacej je przestrzeni.

• Plusy HDBSCAN:
  - dla danych o mocno zroznicowanej gestosci jest lepszy oraz szybszy w porownaniu do DBSCAN
  - HDBSCAN w czasie dzialania odrzuca male odrosty punktow zachowujac najwieksze klastry okreslone przez parametr minimalnego rozmiaru klastra - dzieki temu rowniez dendogram algorytmu HDBSCAN jest mniej skomplikowany
• Minusy HDBSCAN:
  - mimo szybszego czasu dzialania jest duzo ciezszy w zrozumieniu, przez wiele roznych operacji, ktore algorytm wykonuje w czasie dzialania
• Zastosowanie HDBSCAN z parametrem minPts wynoszacym odpowiednio:
10, 50, 100, 200
hdbscan_cluster1 <- hdbscan(data, minPts = 10)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster1$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=10") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

hdbscan_cluster2 <- hdbscan(data, minPts = 50)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster2$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=50") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

hdbscan_cluster3 <- hdbscan(data, minPts = 100)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster3$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=100") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

hdbscan_cluster4 <- hdbscan(data, minPts = 150)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = hdbscan_cluster4$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("HDBSCAN, minPts=150") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

• Wnioski z HDBSCAN:
1) Przy parametrze minPts = 10 dominuje duzy klaster obejmujacy dzielnice Stare Miasto, a takze poludniowe dzielnice sasiadujace takie jak Podgorze, 
Lagiewniki, wschodnie Grzegorzki, Czyzyny, Bienczyce, Mistrzejowice, oraz polnocne Pradnik Bialy i Czerwony.

2) Przy parametrze minPts = 50 widac juz dwa klastry, pierwszy najwiekszy obejmujacy Stare Miasto i okolice, drugi w obszarze Mistrzejowic, Bienczyce oraz Nowa Hute.

3) Dla parametru minPts = 100 brak wiekszych zmian w porownaniu do wyniku z minPts = 50.

4) Dla parametru minPts = 150 brak klastrow.
optics_out1 <- optics(data, minPts = 10)
optics_cluster1 <- extractDBSCAN(optics_out1, eps_cl = 10)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster1$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=10") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

optics_out2 <- optics(data, minPts = 10)
optics_cluster2 <- extractDBSCAN(optics_out2, eps_cl = 50)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster2$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=50") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

optics_out3 <- optics(data, minPts = 10)
optics_cluster3 <- extractDBSCAN(optics_out3, eps_cl = 100)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster3$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=100") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

optics_out4 <- optics(data, minPts = 10)
optics_cluster4 <- extractDBSCAN(optics_out4, eps_cl = 650)
map_cracow + geom_point(data = data, aes(x = Lon, y = Lat, color = optics_cluster4$cluster), size = 1.2) + scale_colour_viridis_c(option = "magma", name = "Intensywnosc") + ggtitle("OPTICS, minPts=10, eps_cl=650") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Dlugosc geograficzna", y = "Szerokosc geograficzna")

• OPTICS (ordering points to identify the clustering structure) - algorytm wykorzystujacy drzewa kd, stworzony w celu wyeliminowania wad algorytmu DBSCAN. OPTICS w trakcie dzialania dynamicznie rozszerza promien wyszukiwania dookola kazdego przypadku, zamiast ustalania z gory okreslonej wartosci promienia. Po zakonczeniu dzialania, a wie po przejsciu wszystkich przypadkow zwracana jest kolejnosc przetwarzania - "odwiedzin" oraz wynik osiagalnosci.

• Zalety:
- brak narzucania z gory wielkosci promienia
- nie wymaga parametrow gestosci
- kolejnosc klastrowania moze byc uzyteczna przy wyodrebnianiu informacji o klastrowaniu

• Wady:
- algorytm nie radzi sobie z wielkowymiarowymi danymi
- tworzy tylko porzadek klastrowy